Import packages

library(sf)
library(mapview)
library(dplyr)
library(tidyr)
library(lubridate)
library(evd)  
library(parallel)  
library(FNN)
library(stringr)

Load map

pathshp <- "~/comp9991/data/geo/DOPGrid.shp"
map <- st_read(pathshp, quiet = TRUE)
print(map)
Simple feature collection with 17450 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 274189.9 ymin: 2901792 xmax: 532797 ymax: 3324561
Projected CRS: NAD83 / UTM zone 17N
First 10 features:
   OBJECTID  PIXEL Shape_Leng Shape_Area                       geometry
1     14202 131641   7697.982    3699323 POLYGON ((411365.9 3234204,...
2     14203 131642   7697.949    3699292 POLYGON ((413224.1 3234190,...
3     14204 131643   7697.920    3699264 POLYGON ((415082.4 3234176,...
4     14205 131644   7697.892    3699237 POLYGON ((416940.6 3234162,...
5     14206 131645   7697.862    3699208 POLYGON ((418798.9 3234149,...
6     14207 131646   7697.834    3699181 POLYGON ((420657.1 3234136,...
7     14208 131647   7697.806    3699155 POLYGON ((422515.3 3234123,...
8     14209 131648   7697.780    3699130 POLYGON ((424373.5 3234111,...
9     14210 131649   7697.756    3699106 POLYGON ((426231.8 3234098,...
10    14211 131650   7697.728    3699080 POLYGON ((428090 3234087, 4...
mapview(map)

Load data

data <- read.csv("~/comp9991/data/PixelRain15min_1995.csv",header = TRUE)

change_time_resolution <- function(data, res = "default"){
  # data: dataframe
  data <- data %>%
  mutate(DateTime = dmy_hms(DateTime),    
         Date = as.Date(DateTime),
         Hour = floor_date(DateTime, "hour"))
  
  if (res == "default"){
    return(data)
  }
  else if (res == "hourly"){
    data_hourly <- data %>% 
      group_by(Hour) %>% 
      summarise(across(starts_with("P"), \(x) sum(x, na.rm = TRUE)))
    
    return(data_hourly)
  }
  else if (res == "daily"){
    data_daily <- data %>%
      group_by(Date) %>%
      summarise(across(starts_with("P"), \(x) sum(x, na.rm = TRUE)))
    return(data_daily)
  }
}

hourly

data_hourly <- change_time_resolution(data,res = "hourly")
class(data_hourly)
[1] "tbl_df"     "tbl"        "data.frame"
all_pixel_IDs <- names(data_hourly)[-1]
all_pixel_IDs <- as.numeric(stringi::stri_extract_first(all_pixel_IDs, regex = "[0-9]+"))
sub_map <- map %>%
  filter(PIXEL %in% all_pixel_IDs)

mapview(sub_map)

cut a square: 14*14

cutid1 <- seq(11951+2, 11951 + 13+2)
cutid2 <- seq(11874+2, 11874 + 13+2)
cutid3 <- seq(11796+2, 11796 + 13+2)
cutid4 <- seq(11716+2, 11716 + 13+2)
cutid5 <- seq(11635+2, 11635 + 13+2)
cutid6 <- seq(11553+2, 11553 + 13+2)
cutid7 <- seq(11470+2, 11470 + 13+2)
cutid8 <- seq(11386+2, 11386 + 13+2)
cutid9 <- seq(11301+2, 11301 + 13+2)
cutid10 <- seq(11215+2, 11215 + 13+2)
cutid11 <- seq(11129+2, 11129 + 13+2)
cutid12 <- seq(11042+2, 11042 + 13+2)
cutid13 <- seq(10954+2, 10954 + 13+2)
cutid14 <- seq(10866+2, 10866 + 13+2)

cut_id <- c(
  cutid1,
  cutid2,
  cutid3,
  cutid4,
  cutid5, 
  cutid6,
  cutid7,
  cutid8,
  cutid9,
  cutid10,
  cutid11,
  cutid12,
  cutid13,
  cutid14
)
square_map <- sub_map %>%
  filter(OBJECTID %in% cut_id)

mapview(square_map)
data_hourly <- data_hourly[rowSums(data_hourly[, -1]) > 0, ]
sub_pixel <- square_map$PIXEL
idx <- all_pixel_IDs %in% sub_pixel
idx <- c(TRUE, idx)
data_hourly_sub <- data_hourly[, idx]
val <- data_hourly_sub[20, 2:ncol(data_hourly_sub)]

val_long <- val %>%
  pivot_longer(
    cols = starts_with("P"),
    names_to = "PIXEL_COL_NAME",
    values_to = "value"  
  ) %>%
  mutate(
    PIXEL = as.integer(str_remove(PIXEL_COL_NAME, pattern = "P"))
  ) %>%
  select(PIXEL, value)

square_map_value <- square_map %>%
  left_join(val_long, by = "PIXEL")

mapview(square_map_value, zcol="value")

filter extremes

rain_hourly_mat <- as.matrix(data_hourly_sub[, -1])
rain_hourly_uniform <- apply(rain_hourly_mat, 2, function(col) {
  ind <- col > 0
  u <- numeric(length(col))
  u[ind] <- rank(col[ind]) / (sum(ind) + 1)
  u
})

eps <- 1e-10
rain_hourly_uniform[rain_hourly_uniform <= 0] <- eps
rain_hourly_uniform[rain_hourly_uniform >= 1] <- 1 - eps

rain_hourly_pareto <- qgpd(rain_hourly_uniform, xi = 1, beta = 1, mu = 1)

L1_risk_hourly <- rowMeans(rain_hourly_pareto, na.rm = TRUE)
u_L1_hourly   <- quantile(L1_risk_hourly, 0.90, na.rm = TRUE)
extreme_L1_hourly <- rain_hourly_pareto[L1_risk_hourly > u_L1_hourly, ]

idx of timestep

idx_ext_hourly <- which(L1_risk_hourly > u_L1_hourly & !is.na(L1_risk_hourly))
times_hourly <- data_hourly_sub[[1]][idx_ext_hourly]

compute the coord

square_map_m <- st_transform(square_map, crs = 3577)
centroids_m <- st_centroid(square_map_m)
Warning: st_centroid assumes attributes are constant over geometries
coords_m <- st_coordinates(centroids_m)    
nn <- get.knn(coords_m, k = 2)$nn.dist[, 2]
cell_len_m <- median(nn)

coord <- coords_m / cell_len_m

bayesian mixture model

build_Gamma <- function(coord, alpha, theta){
  D <- as.matrix(dist(coord))
  G <- (D / theta)^alpha
  diag(G) <- 0 # \gamma(0) = 0
  return(G)
}

prep_br <- function(Gamma) {
  d <- nrow(Gamma)
  out <- vector("list", d)
  for (j in 1:d) {
    idx <- setdiff(seq_len(d), j)
    Sigj <- ( outer(Gamma[idx, j], rep(1, d-1)) +
                outer(rep(1, d-1), Gamma[idx, j]) -
                Gamma[idx, idx] )
    Sigj <- Sigj + .Machine$double.eps * diag(d-1)
    Lj   <- chol(Sigj)
    out[[j]] <- list(idx = idx, L = Lj, Gamma_col = Gamma[, j])
  }
  out
}

intensity_logskew <- function(x,par,alpha.para=TRUE,log=TRUE){
  sigma = par[[1]]
  if(!is.matrix(x)){x <- matrix(x,nrow=1)}
  n = ncol(x)
  if(n==1) return(1/(x^2))
  omega2 = diag(sigma)
  chol.sigma = chol(sigma)
  inv.sigma = chol2inv(chol.sigma)
  logdet.sigma = sum(log(diag(chol.sigma)))*2
  if(alpha.para){
    alpha = par[[2]]
    delta = c(sigma %*% alpha)/sqrt(c(1+alpha %*% sigma %*% alpha))
  }else{
    delta = par[[2]]
    alpha = c(1 - delta %*% inv.sigma %*% delta)^(-1/2) * c(inv.sigma %*% delta)
  }
  a = log(2) + pnorm(delta,log.p=TRUE)
  q = rowSums(inv.sigma)
  sum.q = sum(q);sum.alpha = sum(alpha)
  q.mat = matrix(q,n,n,byrow=TRUE)
  x.log = log(x)
  x.circ = x.log + matrix(a,nrow=nrow(x),ncol=n,byrow=TRUE)
  beta = (1+sum.alpha^2/sum.q)^(-0.5)
  tau.tilde = apply(x.circ,1,function(x.i)  beta * sum((alpha - sum.alpha*q/sum.q) * (x.i + omega2/2))+ beta*sum.alpha/sum.q)
  A = inv.sigma - q %*% t(q)/sum.q
  val = -(n-3)/2 * log(2) - (n-1)/2*log(pi)-1/2*logdet.sigma - 1/2*log(sum.q) - 1/2 * (sum(q*omega2)-1)/sum.q - 1/8*c(omega2 %*% A %*% omega2) 
  val = val - rowSums(x.log) - 1/2 * apply(x.circ,1,function(x.i) c(x.i %*% A %*% x.i) + sum(x.i * (2*q/sum.q + c(A %*% omega2)))) + pnorm(tau.tilde,log.p=TRUE)
  if(log)
    return(val)
  else
    return(exp(val))    
}

dlog_beta  <- function(x, a, b) dbeta(x, a, b, log = TRUE) # mixture weight prior
dlog_gamma <- function(x, a, b) dgamma(x, shape = a, rate = b, log = TRUE)

component_loglik <- function(X, coord, alpha, theta, anchor = 1){
  Gamma <- build_Gamma(coord = coord, alpha = alpha, theta = theta)
  prep  <- prep_br(Gamma)
  idx   <- prep[[anchor]]$idx
  L     <- prep[[anchor]]$L
  Sigma <- t(L) %*% L
  Xsub  <- as.matrix(X[, idx, drop = FALSE])
  intensity_logskew(
    Xsub, par = list(Sigma, delta = rep(0, ncol(Xsub))),
    alpha.para = FALSE, log = TRUE
  )
}

mix_loglik <- function(X, coord, alpha, th1, th2, anchor = 1) {
  if (!is.finite(th1) || th1 <= 0) return(list(l1 = rep(-Inf, nrow(X)), l2 = rep(-Inf, nrow(X))))
  if (!is.finite(th2) || th2 <= 0) return(list(l1 = rep(-Inf, nrow(X)), l2 = rep(-Inf, nrow(X))))
  l1 <- component_loglik(X, coord, alpha = alpha, theta = th1, anchor = anchor)
  l2 <- component_loglik(X, coord, alpha = alpha, theta = th2, anchor = anchor)
  list(l1 = l1, l2 = l2)
}

tau_from_ll <- function(l1, l2, lambda) {
  a <- log(lambda)     + l1
  b <- log(1 - lambda) + l2
  m <- pmax(a, b)
  exp(a - m) / (exp(a - m) + exp(b - m))
}

mcmc_br <- function(
    data, coord, alpha,
    n_iter   = 4000,
    burn     = 1000,
    thin     = 2,
    # priors
    a_lambda = 1, b_lambda = 1,   # for mixture weight lambda
    a_th2 = 1, b_th2 = 1,         # for theta2
    a_del = 1, b_del = 1,         # for delta
    # proposal sds on log-scale
    prop_sd_log_th2  = 0.10,
    prop_sd_log_del  = 0.10,
    # initials
    init_lambda = 0.5,  # mixture weight
    init_th2    = 15,   # theta2
    init_del    = 20,   # delta
    anchor = 1,
    verbose_every = 100
){
  stopifnot(init_lambda > 0, init_lambda < 1, init_th2 > 0, init_del > 0)
  
  keep_idx <- seq.int(burn + 1, n_iter, by = thin)
  n_keep   <- length(keep_idx)
  
  draws <- data.frame(
    iter     = keep_idx,
    lambda   = NA_real_,   # mixture weight
    theta1   = NA_real_,
    theta2   = NA_real_
  )
  z_draws <- matrix(NA_integer_, nrow = n_keep, ncol = nrow(data))
  tau_bar <- rep(0, nrow(data))
  
  # init
  lambda <- init_lambda
  th2    <- init_th2
  del    <- init_del
  th1    <- th2 + del
  
  ll  <- mix_loglik(X = data, coord = coord, alpha = alpha, th1 = th1, th2 = th2, anchor = anchor)
  tau <- tau_from_ll(ll$l1, ll$l2, lambda = lambda)
  z   <- rbinom(n = length(tau), size = 1, prob = tau)
  
  acc_block <- 0L; n_prop <- 0L
  
  # helper: complete-data log-posterior in η-space (η2=log θ2, ηd=log δ)
  logpost_eta <- function(eta2, etad, z_vec) {
    th2p <- exp(eta2)
    delp <- exp(etad)
    th1p <- th2p + delp
    
    llp <- mix_loglik(X = data, coord = coord, alpha = alpha,
                      th1 = th1p, th2 = th2p, anchor = anchor)
    
    lp_ll    <- sum(z_vec * llp$l1 + (1 - z_vec) * llp$l2)
    lp_prior <- dlog_gamma(th2p, a_th2, b_th2) + dlog_gamma(delp, a_del, b_del)
    lp_jac   <- eta2 + etad
    list(value = lp_ll + lp_prior + lp_jac,
         th1 = th1p, th2 = th2p, del = delp, l1_vec = llp$l1, l2_vec = llp$l2)
  }
  
  keep_cursor <- 1L
  for (it in seq_len(n_iter)) {
    
    ## (1) blocked MH for (θ2, δ) on log-scale
    n_prop <- n_prop + 1L
    eta2  <- log(th2);  etad <- log(del)
    eta2p <- rnorm(1, mean = eta2, sd = prop_sd_log_th2)
    etadp <- rnorm(1, mean = etad, sd = prop_sd_log_del)
    
    # current value via same function
    curr <- logpost_eta(eta2, etad, z)
    prop <- logpost_eta(eta2p, etadp, z)
    
    if (log(runif(1)) < (prop$value - curr$value)) {
      th2 <- prop$th2;  del <- prop$del;  th1 <- prop$th1
      ll$l1 <- prop$l1_vec; ll$l2 <- prop$l2_vec
      acc_block <- acc_block + 1L
    }
    # else keep (th1,th2,del,ll) as-is
    
    ## (2) Gibbs for lambda | z (Beta-conjugate)
    n1 <- sum(z); n0 <- length(z) - n1
    lambda <- rbeta(1, a_lambda + n1, b_lambda + n0)
    
    ## (3) Gibbs for z | x, θ, λ
    tau <- tau_from_ll(ll$l1, ll$l2, lambda)
    z   <- rbinom(length(tau), 1, tau)
    
    ## store thinned draws
    if (it %in% keep_idx) {
      tau_bar <- tau_bar + tau / n_keep
      draws$lambda[keep_cursor] <- lambda
      draws$theta1[keep_cursor] <- th1
      draws$theta2[keep_cursor] <- th2
      z_draws[keep_cursor, ]    <- z
      keep_cursor <- keep_cursor + 1L
    }
    
    if (verbose_every > 0 && it %% verbose_every == 0) {
      cat(sprintf("iter %d | lambda=%.3f θ1=%.3f θ2=%.3f | acc(θ-block)=%.3f\n",
                  it, lambda, th1, th2, acc_block / n_prop))
    }
  }
  
  list(
    draws = draws,
    z_draws = z_draws,
    accept_rate_theta_block = acc_block / n_prop,
    tau_bar = tau_bar,
    hard_labels = ifelse(tau_bar > 0.5, 1, 2),
    post_mean = c(lambda = mean(draws$lambda),
                  theta1 = mean(draws$theta1),
                  theta2 = mean(draws$theta2)),
    post_median = c(lambda = median(draws$lambda),
                    theta1 = median(draws$theta1),
                    theta2 = median(draws$theta2))
  )
}

fit the data

alpha <- 1.5

fit1 <- mcmc_br(
  data = extreme_L1_hourly, coord = coord, alpha = alpha,
  n_iter = 10000, burn = 2000, thin = 2,
  # priors
  a_lambda = 1, b_lambda = 1,   # prior for mixture weight lambda
  a_th2 = 1, b_th2 = 1,         # prior for theta2
  a_del = 1, b_del = 1,         # prior for delta
  # initials
  init_th2  = 2,
  init_del  = 4,
  # proposals (log-scale)
  prop_sd_log_th2  = 0.01,
  prop_sd_log_del  = 0.01,
  anchor = 1,
  verbose_every = 50
)
iter 50 | lambda=0.429 θ1=6.866 θ2=2.087 | acc(θ-block)=0.440
iter 100 | lambda=0.419 θ1=6.927 θ2=2.082 | acc(θ-block)=0.330
iter 150 | lambda=0.421 θ1=6.854 θ2=2.088 | acc(θ-block)=0.287
iter 200 | lambda=0.402 θ1=6.889 θ2=2.100 | acc(θ-block)=0.265
iter 250 | lambda=0.379 θ1=6.896 θ2=2.104 | acc(θ-block)=0.276
iter 300 | lambda=0.393 θ1=6.877 θ2=2.098 | acc(θ-block)=0.297
iter 350 | lambda=0.390 θ1=6.863 θ2=2.081 | acc(θ-block)=0.300
iter 400 | lambda=0.422 θ1=6.822 θ2=2.092 | acc(θ-block)=0.285
iter 450 | lambda=0.363 θ1=6.900 θ2=2.092 | acc(θ-block)=0.280
iter 500 | lambda=0.417 θ1=6.800 θ2=2.094 | acc(θ-block)=0.274
iter 550 | lambda=0.428 θ1=6.818 θ2=2.084 | acc(θ-block)=0.291
iter 600 | lambda=0.376 θ1=6.847 θ2=2.092 | acc(θ-block)=0.288
iter 650 | lambda=0.447 θ1=6.826 θ2=2.092 | acc(θ-block)=0.275
iter 700 | lambda=0.397 θ1=6.877 θ2=2.083 | acc(θ-block)=0.283
iter 750 | lambda=0.388 θ1=6.921 θ2=2.087 | acc(θ-block)=0.285
iter 800 | lambda=0.378 θ1=6.824 θ2=2.088 | acc(θ-block)=0.286
iter 850 | lambda=0.411 θ1=6.806 θ2=2.087 | acc(θ-block)=0.287
iter 900 | lambda=0.426 θ1=6.891 θ2=2.106 | acc(θ-block)=0.289
iter 950 | lambda=0.444 θ1=6.879 θ2=2.102 | acc(θ-block)=0.288
iter 1000 | lambda=0.428 θ1=6.795 θ2=2.092 | acc(θ-block)=0.288
iter 1050 | lambda=0.383 θ1=6.841 θ2=2.084 | acc(θ-block)=0.293
iter 1100 | lambda=0.426 θ1=6.824 θ2=2.083 | acc(θ-block)=0.292
iter 1150 | lambda=0.382 θ1=6.842 θ2=2.098 | acc(θ-block)=0.295
iter 1200 | lambda=0.419 θ1=6.878 θ2=2.094 | acc(θ-block)=0.293
iter 1250 | lambda=0.400 θ1=6.866 θ2=2.111 | acc(θ-block)=0.296
iter 1300 | lambda=0.393 θ1=6.793 θ2=2.089 | acc(θ-block)=0.293
iter 1350 | lambda=0.408 θ1=6.819 θ2=2.099 | acc(θ-block)=0.293
iter 1400 | lambda=0.413 θ1=6.850 θ2=2.090 | acc(θ-block)=0.286
iter 1450 | lambda=0.425 θ1=6.883 θ2=2.094 | acc(θ-block)=0.291
iter 1500 | lambda=0.418 θ1=6.905 θ2=2.086 | acc(θ-block)=0.293
iter 1550 | lambda=0.396 θ1=6.833 θ2=2.091 | acc(θ-block)=0.291
iter 1600 | lambda=0.415 θ1=6.790 θ2=2.091 | acc(θ-block)=0.292
iter 1650 | lambda=0.421 θ1=6.896 θ2=2.081 | acc(θ-block)=0.293
iter 1700 | lambda=0.404 θ1=6.905 θ2=2.093 | acc(θ-block)=0.291
iter 1750 | lambda=0.407 θ1=6.872 θ2=2.084 | acc(θ-block)=0.295
iter 1800 | lambda=0.429 θ1=6.812 θ2=2.103 | acc(θ-block)=0.297
iter 1850 | lambda=0.405 θ1=6.776 θ2=2.095 | acc(θ-block)=0.299
iter 1900 | lambda=0.431 θ1=6.843 θ2=2.088 | acc(θ-block)=0.300
iter 1950 | lambda=0.434 θ1=6.839 θ2=2.095 | acc(θ-block)=0.298
iter 2000 | lambda=0.416 θ1=6.843 θ2=2.085 | acc(θ-block)=0.296
iter 2050 | lambda=0.422 θ1=6.875 θ2=2.093 | acc(θ-block)=0.298
iter 2100 | lambda=0.401 θ1=6.827 θ2=2.089 | acc(θ-block)=0.298
iter 2150 | lambda=0.412 θ1=6.935 θ2=2.082 | acc(θ-block)=0.297
iter 2200 | lambda=0.426 θ1=6.834 θ2=2.096 | acc(θ-block)=0.298
iter 2250 | lambda=0.370 θ1=6.881 θ2=2.090 | acc(θ-block)=0.299
iter 2300 | lambda=0.465 θ1=6.816 θ2=2.096 | acc(θ-block)=0.297
iter 2350 | lambda=0.427 θ1=6.869 θ2=2.088 | acc(θ-block)=0.297
iter 2400 | lambda=0.450 θ1=6.830 θ2=2.095 | acc(θ-block)=0.296
iter 2450 | lambda=0.444 θ1=6.794 θ2=2.086 | acc(θ-block)=0.294
iter 2500 | lambda=0.412 θ1=6.773 θ2=2.082 | acc(θ-block)=0.294
iter 2550 | lambda=0.424 θ1=6.824 θ2=2.107 | acc(θ-block)=0.294
iter 2600 | lambda=0.387 θ1=6.920 θ2=2.109 | acc(θ-block)=0.293
iter 2650 | lambda=0.413 θ1=6.863 θ2=2.099 | acc(θ-block)=0.295
iter 2700 | lambda=0.444 θ1=6.824 θ2=2.091 | acc(θ-block)=0.295
iter 2750 | lambda=0.397 θ1=6.922 θ2=2.106 | acc(θ-block)=0.293
iter 2800 | lambda=0.375 θ1=6.832 θ2=2.100 | acc(θ-block)=0.297
iter 2850 | lambda=0.409 θ1=6.829 θ2=2.091 | acc(θ-block)=0.296
iter 2900 | lambda=0.438 θ1=6.822 θ2=2.081 | acc(θ-block)=0.298
iter 2950 | lambda=0.435 θ1=6.937 θ2=2.100 | acc(θ-block)=0.298
iter 3000 | lambda=0.472 θ1=6.796 θ2=2.085 | acc(θ-block)=0.299
iter 3050 | lambda=0.402 θ1=6.873 θ2=2.094 | acc(θ-block)=0.299
iter 3100 | lambda=0.437 θ1=6.870 θ2=2.094 | acc(θ-block)=0.298
iter 3150 | lambda=0.445 θ1=6.828 θ2=2.093 | acc(θ-block)=0.299
iter 3200 | lambda=0.420 θ1=6.910 θ2=2.088 | acc(θ-block)=0.299
iter 3250 | lambda=0.431 θ1=6.842 θ2=2.097 | acc(θ-block)=0.298
iter 3300 | lambda=0.411 θ1=6.885 θ2=2.088 | acc(θ-block)=0.299
iter 3350 | lambda=0.394 θ1=6.894 θ2=2.098 | acc(θ-block)=0.299
iter 3400 | lambda=0.395 θ1=6.814 θ2=2.107 | acc(θ-block)=0.299
iter 3450 | lambda=0.380 θ1=6.784 θ2=2.092 | acc(θ-block)=0.300
iter 3500 | lambda=0.484 θ1=6.853 θ2=2.087 | acc(θ-block)=0.300
iter 3550 | lambda=0.384 θ1=6.856 θ2=2.095 | acc(θ-block)=0.300
iter 3600 | lambda=0.427 θ1=6.850 θ2=2.084 | acc(θ-block)=0.301
iter 3650 | lambda=0.404 θ1=6.902 θ2=2.097 | acc(θ-block)=0.300
iter 3700 | lambda=0.423 θ1=6.773 θ2=2.089 | acc(θ-block)=0.299
iter 3750 | lambda=0.459 θ1=6.865 θ2=2.110 | acc(θ-block)=0.300
iter 3800 | lambda=0.406 θ1=6.808 θ2=2.088 | acc(θ-block)=0.301
iter 3850 | lambda=0.414 θ1=6.800 θ2=2.097 | acc(θ-block)=0.300
iter 3900 | lambda=0.421 θ1=6.834 θ2=2.097 | acc(θ-block)=0.300
iter 3950 | lambda=0.432 θ1=6.845 θ2=2.104 | acc(θ-block)=0.300
iter 4000 | lambda=0.431 θ1=6.840 θ2=2.095 | acc(θ-block)=0.299
iter 4050 | lambda=0.412 θ1=6.800 θ2=2.090 | acc(θ-block)=0.300
iter 4100 | lambda=0.409 θ1=6.854 θ2=2.096 | acc(θ-block)=0.299
iter 4150 | lambda=0.392 θ1=6.789 θ2=2.090 | acc(θ-block)=0.300
iter 4200 | lambda=0.414 θ1=6.844 θ2=2.097 | acc(θ-block)=0.300
iter 4250 | lambda=0.434 θ1=6.853 θ2=2.077 | acc(θ-block)=0.300
iter 4300 | lambda=0.387 θ1=6.849 θ2=2.095 | acc(θ-block)=0.299
iter 4350 | lambda=0.414 θ1=6.901 θ2=2.097 | acc(θ-block)=0.297
iter 4400 | lambda=0.409 θ1=6.923 θ2=2.086 | acc(θ-block)=0.297
iter 4450 | lambda=0.413 θ1=6.781 θ2=2.091 | acc(θ-block)=0.297
iter 4500 | lambda=0.420 θ1=6.883 θ2=2.095 | acc(θ-block)=0.298
iter 4550 | lambda=0.437 θ1=6.876 θ2=2.091 | acc(θ-block)=0.299
iter 4600 | lambda=0.375 θ1=6.823 θ2=2.099 | acc(θ-block)=0.298
iter 4650 | lambda=0.438 θ1=6.866 θ2=2.091 | acc(θ-block)=0.298
iter 4700 | lambda=0.420 θ1=6.817 θ2=2.094 | acc(θ-block)=0.299
iter 4750 | lambda=0.412 θ1=6.921 θ2=2.094 | acc(θ-block)=0.299
iter 4800 | lambda=0.391 θ1=6.841 θ2=2.093 | acc(θ-block)=0.299
iter 4850 | lambda=0.402 θ1=6.868 θ2=2.100 | acc(θ-block)=0.299
iter 4900 | lambda=0.390 θ1=6.861 θ2=2.100 | acc(θ-block)=0.300
iter 4950 | lambda=0.407 θ1=6.825 θ2=2.090 | acc(θ-block)=0.300
iter 5000 | lambda=0.395 θ1=6.845 θ2=2.096 | acc(θ-block)=0.299
iter 5050 | lambda=0.407 θ1=6.833 θ2=2.090 | acc(θ-block)=0.298
iter 5100 | lambda=0.396 θ1=6.895 θ2=2.091 | acc(θ-block)=0.298
iter 5150 | lambda=0.368 θ1=6.906 θ2=2.109 | acc(θ-block)=0.297
iter 5200 | lambda=0.413 θ1=6.820 θ2=2.086 | acc(θ-block)=0.299
iter 5250 | lambda=0.426 θ1=6.832 θ2=2.099 | acc(θ-block)=0.298
iter 5300 | lambda=0.457 θ1=7.003 θ2=2.097 | acc(θ-block)=0.298
iter 5350 | lambda=0.407 θ1=6.853 θ2=2.087 | acc(θ-block)=0.298
iter 5400 | lambda=0.396 θ1=6.839 θ2=2.095 | acc(θ-block)=0.298
iter 5450 | lambda=0.373 θ1=6.909 θ2=2.092 | acc(θ-block)=0.298
iter 5500 | lambda=0.405 θ1=6.892 θ2=2.080 | acc(θ-block)=0.298
iter 5550 | lambda=0.404 θ1=6.831 θ2=2.089 | acc(θ-block)=0.298
iter 5600 | lambda=0.442 θ1=6.887 θ2=2.088 | acc(θ-block)=0.297
iter 5650 | lambda=0.386 θ1=6.843 θ2=2.078 | acc(θ-block)=0.297
iter 5700 | lambda=0.393 θ1=6.896 θ2=2.096 | acc(θ-block)=0.297
iter 5750 | lambda=0.396 θ1=6.903 θ2=2.101 | acc(θ-block)=0.297
iter 5800 | lambda=0.429 θ1=6.786 θ2=2.094 | acc(θ-block)=0.297
iter 5850 | lambda=0.437 θ1=6.881 θ2=2.083 | acc(θ-block)=0.297
iter 5900 | lambda=0.443 θ1=6.841 θ2=2.098 | acc(θ-block)=0.297
iter 5950 | lambda=0.409 θ1=6.806 θ2=2.096 | acc(θ-block)=0.297
iter 6000 | lambda=0.426 θ1=6.859 θ2=2.097 | acc(θ-block)=0.297
iter 6050 | lambda=0.379 θ1=6.868 θ2=2.092 | acc(θ-block)=0.298
iter 6100 | lambda=0.402 θ1=6.841 θ2=2.103 | acc(θ-block)=0.298
iter 6150 | lambda=0.427 θ1=6.837 θ2=2.093 | acc(θ-block)=0.299
iter 6200 | lambda=0.413 θ1=6.866 θ2=2.069 | acc(θ-block)=0.300
iter 6250 | lambda=0.425 θ1=6.815 θ2=2.083 | acc(θ-block)=0.300
iter 6300 | lambda=0.459 θ1=6.787 θ2=2.097 | acc(θ-block)=0.300
iter 6350 | lambda=0.415 θ1=6.884 θ2=2.107 | acc(θ-block)=0.300
iter 6400 | lambda=0.355 θ1=6.890 θ2=2.096 | acc(θ-block)=0.300
iter 6450 | lambda=0.413 θ1=6.852 θ2=2.079 | acc(θ-block)=0.300
iter 6500 | lambda=0.429 θ1=6.851 θ2=2.105 | acc(θ-block)=0.301
iter 6550 | lambda=0.449 θ1=6.868 θ2=2.075 | acc(θ-block)=0.300
iter 6600 | lambda=0.399 θ1=6.870 θ2=2.098 | acc(θ-block)=0.300
iter 6650 | lambda=0.432 θ1=6.878 θ2=2.107 | acc(θ-block)=0.300
iter 6700 | lambda=0.423 θ1=6.850 θ2=2.099 | acc(θ-block)=0.300
iter 6750 | lambda=0.437 θ1=6.833 θ2=2.099 | acc(θ-block)=0.299
iter 6800 | lambda=0.386 θ1=6.915 θ2=2.101 | acc(θ-block)=0.300
iter 6850 | lambda=0.443 θ1=6.844 θ2=2.093 | acc(θ-block)=0.300
iter 6900 | lambda=0.398 θ1=6.832 θ2=2.103 | acc(θ-block)=0.299
iter 6950 | lambda=0.385 θ1=6.863 θ2=2.100 | acc(θ-block)=0.299
iter 7000 | lambda=0.353 θ1=6.860 θ2=2.104 | acc(θ-block)=0.299
iter 7050 | lambda=0.455 θ1=6.841 θ2=2.089 | acc(θ-block)=0.299
iter 7100 | lambda=0.396 θ1=6.833 θ2=2.095 | acc(θ-block)=0.298
iter 7150 | lambda=0.457 θ1=6.860 θ2=2.080 | acc(θ-block)=0.299
iter 7200 | lambda=0.412 θ1=6.890 θ2=2.101 | acc(θ-block)=0.299
iter 7250 | lambda=0.403 θ1=6.835 θ2=2.089 | acc(θ-block)=0.299
iter 7300 | lambda=0.458 θ1=6.819 θ2=2.083 | acc(θ-block)=0.299
iter 7350 | lambda=0.385 θ1=6.919 θ2=2.083 | acc(θ-block)=0.300
iter 7400 | lambda=0.426 θ1=6.829 θ2=2.109 | acc(θ-block)=0.299
iter 7450 | lambda=0.454 θ1=6.874 θ2=2.100 | acc(θ-block)=0.299
iter 7500 | lambda=0.460 θ1=6.803 θ2=2.097 | acc(θ-block)=0.299
iter 7550 | lambda=0.413 θ1=6.808 θ2=2.104 | acc(θ-block)=0.298
iter 7600 | lambda=0.453 θ1=6.820 θ2=2.092 | acc(θ-block)=0.298
iter 7650 | lambda=0.386 θ1=6.839 θ2=2.098 | acc(θ-block)=0.298
iter 7700 | lambda=0.418 θ1=6.844 θ2=2.098 | acc(θ-block)=0.299
iter 7750 | lambda=0.376 θ1=6.866 θ2=2.082 | acc(θ-block)=0.300
iter 7800 | lambda=0.437 θ1=6.842 θ2=2.083 | acc(θ-block)=0.300
iter 7850 | lambda=0.366 θ1=6.843 θ2=2.085 | acc(θ-block)=0.299
iter 7900 | lambda=0.479 θ1=6.838 θ2=2.094 | acc(θ-block)=0.299
iter 7950 | lambda=0.417 θ1=6.795 θ2=2.093 | acc(θ-block)=0.300
iter 8000 | lambda=0.437 θ1=6.837 θ2=2.081 | acc(θ-block)=0.300
iter 8050 | lambda=0.422 θ1=6.945 θ2=2.096 | acc(θ-block)=0.300
iter 8100 | lambda=0.424 θ1=6.814 θ2=2.091 | acc(θ-block)=0.299
iter 8150 | lambda=0.451 θ1=6.873 θ2=2.100 | acc(θ-block)=0.299
iter 8200 | lambda=0.439 θ1=6.843 θ2=2.085 | acc(θ-block)=0.299
iter 8250 | lambda=0.427 θ1=6.859 θ2=2.099 | acc(θ-block)=0.299
iter 8300 | lambda=0.430 θ1=6.893 θ2=2.094 | acc(θ-block)=0.298
iter 8350 | lambda=0.424 θ1=6.869 θ2=2.102 | acc(θ-block)=0.298
iter 8400 | lambda=0.390 θ1=6.865 θ2=2.085 | acc(θ-block)=0.298
iter 8450 | lambda=0.393 θ1=6.830 θ2=2.092 | acc(θ-block)=0.299
iter 8500 | lambda=0.430 θ1=6.809 θ2=2.077 | acc(θ-block)=0.299
iter 8550 | lambda=0.416 θ1=6.837 θ2=2.100 | acc(θ-block)=0.299
iter 8600 | lambda=0.402 θ1=6.826 θ2=2.085 | acc(θ-block)=0.300
iter 8650 | lambda=0.465 θ1=6.829 θ2=2.091 | acc(θ-block)=0.300
iter 8700 | lambda=0.415 θ1=6.816 θ2=2.096 | acc(θ-block)=0.300
iter 8750 | lambda=0.392 θ1=6.839 θ2=2.081 | acc(θ-block)=0.300
iter 8800 | lambda=0.433 θ1=6.839 θ2=2.099 | acc(θ-block)=0.299
iter 8850 | lambda=0.426 θ1=6.795 θ2=2.091 | acc(θ-block)=0.298
iter 8900 | lambda=0.430 θ1=6.912 θ2=2.083 | acc(θ-block)=0.298
iter 8950 | lambda=0.391 θ1=6.844 θ2=2.099 | acc(θ-block)=0.298
iter 9000 | lambda=0.433 θ1=6.830 θ2=2.102 | acc(θ-block)=0.298
iter 9050 | lambda=0.420 θ1=6.875 θ2=2.101 | acc(θ-block)=0.298
iter 9100 | lambda=0.401 θ1=6.876 θ2=2.093 | acc(θ-block)=0.297
iter 9150 | lambda=0.402 θ1=6.886 θ2=2.107 | acc(θ-block)=0.297
iter 9200 | lambda=0.428 θ1=6.846 θ2=2.092 | acc(θ-block)=0.296
iter 9250 | lambda=0.401 θ1=6.897 θ2=2.094 | acc(θ-block)=0.295
iter 9300 | lambda=0.413 θ1=6.791 θ2=2.093 | acc(θ-block)=0.295
iter 9350 | lambda=0.377 θ1=6.864 θ2=2.098 | acc(θ-block)=0.295
iter 9400 | lambda=0.418 θ1=6.830 θ2=2.100 | acc(θ-block)=0.296
iter 9450 | lambda=0.397 θ1=6.901 θ2=2.081 | acc(θ-block)=0.296
iter 9500 | lambda=0.363 θ1=6.867 θ2=2.092 | acc(θ-block)=0.295
iter 9550 | lambda=0.384 θ1=6.833 θ2=2.092 | acc(θ-block)=0.295
iter 9600 | lambda=0.405 θ1=6.873 θ2=2.087 | acc(θ-block)=0.295
iter 9650 | lambda=0.390 θ1=6.851 θ2=2.099 | acc(θ-block)=0.296
iter 9700 | lambda=0.423 θ1=6.871 θ2=2.081 | acc(θ-block)=0.296
iter 9750 | lambda=0.448 θ1=6.902 θ2=2.092 | acc(θ-block)=0.296
iter 9800 | lambda=0.417 θ1=6.867 θ2=2.097 | acc(θ-block)=0.296
iter 9850 | lambda=0.419 θ1=6.903 θ2=2.093 | acc(θ-block)=0.295
iter 9900 | lambda=0.435 θ1=6.832 θ2=2.102 | acc(θ-block)=0.295
iter 9950 | lambda=0.488 θ1=6.849 θ2=2.091 | acc(θ-block)=0.295
iter 10000 | lambda=0.423 θ1=6.845 θ2=2.100 | acc(θ-block)=0.295
save(fit1, file = "application_fit_hourly_final.RData")
df <- fit1$draws
g1 <- ggplot(df, aes(iter, theta1)) + geom_line() +
  labs(title = expression(trace~theta[1]))
g2 <- ggplot(df, aes(iter, theta2)) + geom_line() +
  labs(title = expression(trace~theta[2]))
g3 <- ggplot(df, aes(iter, lambda)) + geom_line() +
  labs(title = expression(trace~lambda))
gd1 <- ggplot(df, aes(theta1)) + geom_density() +
  labs(title = expression(density~theta[1]))
gd2 <- ggplot(df, aes(theta2)) + geom_density() +
  labs(title = expression(density~theta[2]))
gd3 <- ggplot(df, aes(lambda)) + geom_density() +
  labs(title = expression(density~lambda))

gridExtra::grid.arrange(g1,g2,g3, gd1,gd2,gd3, ncol = 3)

fit1$post_mean
   lambda    theta1    theta2 
0.4186907 6.8560784 2.0937366 
sd(fit1$draws$theta1)
[1] 0.03948329
sd(fit1$draws$theta2)
[1] 0.008061415
plot_extreme_on_map <- function(i, extreme_mat = extreme_L1_hourly,
                                map_sf = square_map, data_wide = data_hourly_sub) {
  stopifnot(is.matrix(extreme_mat),
            i >= 1, i <= nrow(extreme_mat))
  
  # column order of pixels in the wide daily data used to build rain_mat/extreme_L1
  pix_order <- as.integer(sub("P", "", colnames(data_wide)[-1]))
  
  # values for the i-th extreme event
  vals_i <- as.numeric(extreme_mat[i, ])
  df_i <- data.frame(PIXEL = pix_order, val_i = vals_i)
  
  # join to geometry
  joined <- map_sf %>%
    left_join(df_i, by = "PIXEL")
  
  # quick stretch for nicer legend breaks
  brks <- quantile(df_i$val_i, probs = seq(0, 1, length.out = 6), na.rm = TRUE)
  brks_fixed <- sort(unique(brks))
  # plot
  mapview(
    joined,
    zcol       = "val_i",
    layer.name = paste0("Extreme L1 snapshot (row ", i, ")"),
    at         = brks_fixed
  )
}
plot_extreme_on_map(1)
plot_extreme_on_map(5)

\(\chi^2(S_A, S_B)\) test

chi_q <- function(X, q = 0.98) {
  n <- ncol(X)
  m <- nrow(X)
  
  # thresholds at site level
  u <- apply(X, 2, quantile, probs = q, na.rm = TRUE)
  
  # initialize matrix
  chi_hat <- matrix(NA, n, n)
  
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      Ai <- X[, i] > u[i]
      Aj <- X[, j] > u[j]
      chi_hat[i, j] <- mean(Ai & Aj) / mean(Ai)
      chi_hat[j, i] <- chi_hat[i, j]
    }
  }
  diag(chi_hat) <- 1
  chi_hat
}

dist_matrix <- as.matrix(dist(coord))

stopifnot(length(fit1$hard_labels) == nrow(extreme_L1_hourly))
labs <- fit1$hard_labels

X1 <- extreme_L1_hourly[labs == 1, , drop = FALSE]
X2 <- extreme_L1_hourly[labs == 2, , drop = FALSE]

# 4) compute χ̂_q for both comps and two q levels
chi_95_c1 <- chi_q(X1, q = 0.95)
chi_95_c2 <- chi_q(X2, q = 0.95)
chi_90_c1 <- chi_q(X1, q = 0.90)
chi_90_c2 <- chi_q(X2, q = 0.90)
upper_to_df <- function(chi_mat, dist_mat, q_label, comp_label) {
  idx <- which(upper.tri(chi_mat), arr.ind = TRUE)
  tibble(
    dist = dist_mat[idx],
    chi  = chi_mat[idx],
    q    = q_label,
    comp = comp_label
  ) %>% filter(is.finite(chi))
}

df_all <- bind_rows(
  upper_to_df(chi_95_c1, dist_matrix, "q=0.95", "Comp1"),
  upper_to_df(chi_95_c2, dist_matrix, "q=0.95", "Comp2"),
  upper_to_df(chi_90_c1, dist_matrix, "q=0.90", "Comp1"),
  upper_to_df(chi_90_c2, dist_matrix, "q=0.90", "Comp2"),

)

# 6) plot (mirrors your original style)
ggplot(df_all, aes(x = dist, y = chi, color = comp)) +
  geom_point(alpha = 0.4, size = 0.6) +
  facet_wrap(~ q, scales = "free_y") +
  theme_bw() +
  labs(
    x = expression(h(s[A], s[B])),
    y = expression(hat(chi)[q](s[A], s[B])),
    title = "Empirical " %+% expression(hat(chi)[q]) %+% " vs Distance by MCMC Hard-Label Component"
  ) +
  scale_color_manual(values = c("Comp1" = "blue", "Comp2" = "red"))
for (i in seq_len(ncol(fit1$z_draws))) {
  plot(
    fit1$z_draws[, i],
    type = "l",
    xlab = "iterations",
    ylab = "z draw",
    main = paste("sample", i) # <--- Modified line
  )
}



plot(times_hourly, fit1$tau_bar, type = "l",
     xlab = "Time", ylab = "Posterior mean of z",
     main = "Posterior mean with season")

LS0tCnRpdGxlOiAiQXBwbGljYXRpb24gRmluYWwgVmVyc2lvbiIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyBJbXBvcnQgcGFja2FnZXMKYGBge3J9CmxpYnJhcnkoc2YpCmxpYnJhcnkobWFwdmlldykKbGlicmFyeShkcGx5cikKbGlicmFyeSh0aWR5cikKbGlicmFyeShsdWJyaWRhdGUpCmxpYnJhcnkoZXZkKSAgCmxpYnJhcnkocGFyYWxsZWwpICAKbGlicmFyeShGTk4pCmxpYnJhcnkoc3RyaW5ncikKYGBgCgojIExvYWQgbWFwCmBgYHtyfQpwYXRoc2hwIDwtICJ+L2NvbXA5OTkxL2RhdGEvZ2VvL0RPUEdyaWQuc2hwIgptYXAgPC0gc3RfcmVhZChwYXRoc2hwLCBxdWlldCA9IFRSVUUpCnByaW50KG1hcCkKbWFwdmlldyhtYXApCmBgYAojIExvYWQgZGF0YQpgYGB7cn0KZGF0YSA8LSByZWFkLmNzdigifi9jb21wOTk5MS9kYXRhL1BpeGVsUmFpbjE1bWluXzE5OTUuY3N2IixoZWFkZXIgPSBUUlVFKQoKY2hhbmdlX3RpbWVfcmVzb2x1dGlvbiA8LSBmdW5jdGlvbihkYXRhLCByZXMgPSAiZGVmYXVsdCIpewogICMgZGF0YTogZGF0YWZyYW1lCiAgZGF0YSA8LSBkYXRhICU+JQogIG11dGF0ZShEYXRlVGltZSA9IGRteV9obXMoRGF0ZVRpbWUpLCAgICAKICAgICAgICAgRGF0ZSA9IGFzLkRhdGUoRGF0ZVRpbWUpLAogICAgICAgICBIb3VyID0gZmxvb3JfZGF0ZShEYXRlVGltZSwgImhvdXIiKSkKICAKICBpZiAocmVzID09ICJkZWZhdWx0Iil7CiAgICByZXR1cm4oZGF0YSkKICB9CiAgZWxzZSBpZiAocmVzID09ICJob3VybHkiKXsKICAgIGRhdGFfaG91cmx5IDwtIGRhdGEgJT4lIAogICAgICBncm91cF9ieShIb3VyKSAlPiUgCiAgICAgIHN1bW1hcmlzZShhY3Jvc3Moc3RhcnRzX3dpdGgoIlAiKSwgXCh4KSBzdW0oeCwgbmEucm0gPSBUUlVFKSkpCiAgICAKICAgIHJldHVybihkYXRhX2hvdXJseSkKICB9CiAgZWxzZSBpZiAocmVzID09ICJkYWlseSIpewogICAgZGF0YV9kYWlseSA8LSBkYXRhICU+JQogICAgICBncm91cF9ieShEYXRlKSAlPiUKICAgICAgc3VtbWFyaXNlKGFjcm9zcyhzdGFydHNfd2l0aCgiUCIpLCBcKHgpIHN1bSh4LCBuYS5ybSA9IFRSVUUpKSkKICAgIHJldHVybihkYXRhX2RhaWx5KQogIH0KfQoKCmBgYAoKaG91cmx5CmBgYHtyfQpkYXRhX2hvdXJseSA8LSBjaGFuZ2VfdGltZV9yZXNvbHV0aW9uKGRhdGEscmVzID0gImhvdXJseSIpCmNsYXNzKGRhdGFfaG91cmx5KQpgYGAKCmBgYHtyfQphbGxfcGl4ZWxfSURzIDwtIG5hbWVzKGRhdGFfaG91cmx5KVstMV0KYWxsX3BpeGVsX0lEcyA8LSBhcy5udW1lcmljKHN0cmluZ2k6OnN0cmlfZXh0cmFjdF9maXJzdChhbGxfcGl4ZWxfSURzLCByZWdleCA9ICJbMC05XSsiKSkKYGBgCgpgYGB7cn0Kc3ViX21hcCA8LSBtYXAgJT4lCiAgZmlsdGVyKFBJWEVMICVpbiUgYWxsX3BpeGVsX0lEcykKCm1hcHZpZXcoc3ViX21hcCkKYGBgCgpjdXQgYSBzcXVhcmU6IDE0KjE0CmBgYHtyfQpjdXRpZDEgPC0gc2VxKDExOTUxKzIsIDExOTUxICsgMTMrMikKY3V0aWQyIDwtIHNlcSgxMTg3NCsyLCAxMTg3NCArIDEzKzIpCmN1dGlkMyA8LSBzZXEoMTE3OTYrMiwgMTE3OTYgKyAxMysyKQpjdXRpZDQgPC0gc2VxKDExNzE2KzIsIDExNzE2ICsgMTMrMikKY3V0aWQ1IDwtIHNlcSgxMTYzNSsyLCAxMTYzNSArIDEzKzIpCmN1dGlkNiA8LSBzZXEoMTE1NTMrMiwgMTE1NTMgKyAxMysyKQpjdXRpZDcgPC0gc2VxKDExNDcwKzIsIDExNDcwICsgMTMrMikKY3V0aWQ4IDwtIHNlcSgxMTM4NisyLCAxMTM4NiArIDEzKzIpCmN1dGlkOSA8LSBzZXEoMTEzMDErMiwgMTEzMDEgKyAxMysyKQpjdXRpZDEwIDwtIHNlcSgxMTIxNSsyLCAxMTIxNSArIDEzKzIpCmN1dGlkMTEgPC0gc2VxKDExMTI5KzIsIDExMTI5ICsgMTMrMikKY3V0aWQxMiA8LSBzZXEoMTEwNDIrMiwgMTEwNDIgKyAxMysyKQpjdXRpZDEzIDwtIHNlcSgxMDk1NCsyLCAxMDk1NCArIDEzKzIpCmN1dGlkMTQgPC0gc2VxKDEwODY2KzIsIDEwODY2ICsgMTMrMikKCmN1dF9pZCA8LSBjKAogIGN1dGlkMSwKICBjdXRpZDIsCiAgY3V0aWQzLAogIGN1dGlkNCwKICBjdXRpZDUsIAogIGN1dGlkNiwKICBjdXRpZDcsCiAgY3V0aWQ4LAogIGN1dGlkOSwKICBjdXRpZDEwLAogIGN1dGlkMTEsCiAgY3V0aWQxMiwKICBjdXRpZDEzLAogIGN1dGlkMTQKKQpgYGAKCmBgYHtyfQpzcXVhcmVfbWFwIDwtIHN1Yl9tYXAgJT4lCiAgZmlsdGVyKE9CSkVDVElEICVpbiUgY3V0X2lkKQoKbWFwdmlldyhzcXVhcmVfbWFwKQpgYGAKCmBgYHtyfQpkYXRhX2hvdXJseSA8LSBkYXRhX2hvdXJseVtyb3dTdW1zKGRhdGFfaG91cmx5WywgLTFdKSA+IDAsIF0Kc3ViX3BpeGVsIDwtIHNxdWFyZV9tYXAkUElYRUwKaWR4IDwtIGFsbF9waXhlbF9JRHMgJWluJSBzdWJfcGl4ZWwKaWR4IDwtIGMoVFJVRSwgaWR4KQpkYXRhX2hvdXJseV9zdWIgPC0gZGF0YV9ob3VybHlbLCBpZHhdCmBgYAoKYGBge3J9CnZhbCA8LSBkYXRhX2hvdXJseV9zdWJbMjAsIDI6bmNvbChkYXRhX2hvdXJseV9zdWIpXQoKdmFsX2xvbmcgPC0gdmFsICU+JQogIHBpdm90X2xvbmdlcigKICAgIGNvbHMgPSBzdGFydHNfd2l0aCgiUCIpLAogICAgbmFtZXNfdG8gPSAiUElYRUxfQ09MX05BTUUiLAogICAgdmFsdWVzX3RvID0gInZhbHVlIiAgCiAgKSAlPiUKICBtdXRhdGUoCiAgICBQSVhFTCA9IGFzLmludGVnZXIoc3RyX3JlbW92ZShQSVhFTF9DT0xfTkFNRSwgcGF0dGVybiA9ICJQIikpCiAgKSAlPiUKICBzZWxlY3QoUElYRUwsIHZhbHVlKQoKc3F1YXJlX21hcF92YWx1ZSA8LSBzcXVhcmVfbWFwICU+JQogIGxlZnRfam9pbih2YWxfbG9uZywgYnkgPSAiUElYRUwiKQoKbWFwdmlldyhzcXVhcmVfbWFwX3ZhbHVlLCB6Y29sPSJ2YWx1ZSIpCmBgYAoKZmlsdGVyIGV4dHJlbWVzCmBgYHtyfQpyYWluX2hvdXJseV9tYXQgPC0gYXMubWF0cml4KGRhdGFfaG91cmx5X3N1YlssIC0xXSkKcmFpbl9ob3VybHlfdW5pZm9ybSA8LSBhcHBseShyYWluX2hvdXJseV9tYXQsIDIsIGZ1bmN0aW9uKGNvbCkgewogIGluZCA8LSBjb2wgPiAwCiAgdSA8LSBudW1lcmljKGxlbmd0aChjb2wpKQogIHVbaW5kXSA8LSByYW5rKGNvbFtpbmRdKSAvIChzdW0oaW5kKSArIDEpCiAgdQp9KQoKZXBzIDwtIDFlLTEwCnJhaW5faG91cmx5X3VuaWZvcm1bcmFpbl9ob3VybHlfdW5pZm9ybSA8PSAwXSA8LSBlcHMKcmFpbl9ob3VybHlfdW5pZm9ybVtyYWluX2hvdXJseV91bmlmb3JtID49IDFdIDwtIDEgLSBlcHMKCnJhaW5faG91cmx5X3BhcmV0byA8LSBxZ3BkKHJhaW5faG91cmx5X3VuaWZvcm0sIHhpID0gMSwgYmV0YSA9IDEsIG11ID0gMSkKCkwxX3Jpc2tfaG91cmx5IDwtIHJvd01lYW5zKHJhaW5faG91cmx5X3BhcmV0bywgbmEucm0gPSBUUlVFKQp1X0wxX2hvdXJseSAgIDwtIHF1YW50aWxlKEwxX3Jpc2tfaG91cmx5LCAwLjkwLCBuYS5ybSA9IFRSVUUpCmV4dHJlbWVfTDFfaG91cmx5IDwtIHJhaW5faG91cmx5X3BhcmV0b1tMMV9yaXNrX2hvdXJseSA+IHVfTDFfaG91cmx5LCBdCmBgYAoKaWR4IG9mIHRpbWVzdGVwCmBgYHtyfQppZHhfZXh0X2hvdXJseSA8LSB3aGljaChMMV9yaXNrX2hvdXJseSA+IHVfTDFfaG91cmx5ICYgIWlzLm5hKEwxX3Jpc2tfaG91cmx5KSkKdGltZXNfaG91cmx5IDwtIGRhdGFfaG91cmx5X3N1YltbMV1dW2lkeF9leHRfaG91cmx5XQpgYGAKCiMgY29tcHV0ZSB0aGUgY29vcmQKYGBge3J9CnNxdWFyZV9tYXBfbSA8LSBzdF90cmFuc2Zvcm0oc3F1YXJlX21hcCwgY3JzID0gMzU3NykKY2VudHJvaWRzX20gPC0gc3RfY2VudHJvaWQoc3F1YXJlX21hcF9tKQpjb29yZHNfbSA8LSBzdF9jb29yZGluYXRlcyhjZW50cm9pZHNfbSkgICAgCmBgYAoKCmBgYHtyfQpubiA8LSBnZXQua25uKGNvb3Jkc19tLCBrID0gMikkbm4uZGlzdFssIDJdCmNlbGxfbGVuX20gPC0gbWVkaWFuKG5uKQoKY29vcmQgPC0gY29vcmRzX20gLyBjZWxsX2xlbl9tCmBgYAoKIyBiYXllc2lhbiBtaXh0dXJlIG1vZGVsCmBgYHtyfQpidWlsZF9HYW1tYSA8LSBmdW5jdGlvbihjb29yZCwgYWxwaGEsIHRoZXRhKXsKICBEIDwtIGFzLm1hdHJpeChkaXN0KGNvb3JkKSkKICBHIDwtIChEIC8gdGhldGEpXmFscGhhCiAgZGlhZyhHKSA8LSAwICMgXGdhbW1hKDApID0gMAogIHJldHVybihHKQp9CgpwcmVwX2JyIDwtIGZ1bmN0aW9uKEdhbW1hKSB7CiAgZCA8LSBucm93KEdhbW1hKQogIG91dCA8LSB2ZWN0b3IoImxpc3QiLCBkKQogIGZvciAoaiBpbiAxOmQpIHsKICAgIGlkeCA8LSBzZXRkaWZmKHNlcV9sZW4oZCksIGopCiAgICBTaWdqIDwtICggb3V0ZXIoR2FtbWFbaWR4LCBqXSwgcmVwKDEsIGQtMSkpICsKICAgICAgICAgICAgICAgIG91dGVyKHJlcCgxLCBkLTEpLCBHYW1tYVtpZHgsIGpdKSAtCiAgICAgICAgICAgICAgICBHYW1tYVtpZHgsIGlkeF0gKQogICAgU2lnaiA8LSBTaWdqICsgLk1hY2hpbmUkZG91YmxlLmVwcyAqIGRpYWcoZC0xKQogICAgTGogICA8LSBjaG9sKFNpZ2opCiAgICBvdXRbW2pdXSA8LSBsaXN0KGlkeCA9IGlkeCwgTCA9IExqLCBHYW1tYV9jb2wgPSBHYW1tYVssIGpdKQogIH0KICBvdXQKfQoKaW50ZW5zaXR5X2xvZ3NrZXcgPC0gZnVuY3Rpb24oeCxwYXIsYWxwaGEucGFyYT1UUlVFLGxvZz1UUlVFKXsKICBzaWdtYSA9IHBhcltbMV1dCiAgaWYoIWlzLm1hdHJpeCh4KSl7eCA8LSBtYXRyaXgoeCxucm93PTEpfQogIG4gPSBuY29sKHgpCiAgaWYobj09MSkgcmV0dXJuKDEvKHheMikpCiAgb21lZ2EyID0gZGlhZyhzaWdtYSkKICBjaG9sLnNpZ21hID0gY2hvbChzaWdtYSkKICBpbnYuc2lnbWEgPSBjaG9sMmludihjaG9sLnNpZ21hKQogIGxvZ2RldC5zaWdtYSA9IHN1bShsb2coZGlhZyhjaG9sLnNpZ21hKSkpKjIKICBpZihhbHBoYS5wYXJhKXsKICAgIGFscGhhID0gcGFyW1syXV0KICAgIGRlbHRhID0gYyhzaWdtYSAlKiUgYWxwaGEpL3NxcnQoYygxK2FscGhhICUqJSBzaWdtYSAlKiUgYWxwaGEpKQogIH1lbHNlewogICAgZGVsdGEgPSBwYXJbWzJdXQogICAgYWxwaGEgPSBjKDEgLSBkZWx0YSAlKiUgaW52LnNpZ21hICUqJSBkZWx0YSleKC0xLzIpICogYyhpbnYuc2lnbWEgJSolIGRlbHRhKQogIH0KICBhID0gbG9nKDIpICsgcG5vcm0oZGVsdGEsbG9nLnA9VFJVRSkKICBxID0gcm93U3VtcyhpbnYuc2lnbWEpCiAgc3VtLnEgPSBzdW0ocSk7c3VtLmFscGhhID0gc3VtKGFscGhhKQogIHEubWF0ID0gbWF0cml4KHEsbixuLGJ5cm93PVRSVUUpCiAgeC5sb2cgPSBsb2coeCkKICB4LmNpcmMgPSB4LmxvZyArIG1hdHJpeChhLG5yb3c9bnJvdyh4KSxuY29sPW4sYnlyb3c9VFJVRSkKICBiZXRhID0gKDErc3VtLmFscGhhXjIvc3VtLnEpXigtMC41KQogIHRhdS50aWxkZSA9IGFwcGx5KHguY2lyYywxLGZ1bmN0aW9uKHguaSkgIGJldGEgKiBzdW0oKGFscGhhIC0gc3VtLmFscGhhKnEvc3VtLnEpICogKHguaSArIG9tZWdhMi8yKSkrIGJldGEqc3VtLmFscGhhL3N1bS5xKQogIEEgPSBpbnYuc2lnbWEgLSBxICUqJSB0KHEpL3N1bS5xCiAgdmFsID0gLShuLTMpLzIgKiBsb2coMikgLSAobi0xKS8yKmxvZyhwaSktMS8yKmxvZ2RldC5zaWdtYSAtIDEvMipsb2coc3VtLnEpIC0gMS8yICogKHN1bShxKm9tZWdhMiktMSkvc3VtLnEgLSAxLzgqYyhvbWVnYTIgJSolIEEgJSolIG9tZWdhMikgCiAgdmFsID0gdmFsIC0gcm93U3Vtcyh4LmxvZykgLSAxLzIgKiBhcHBseSh4LmNpcmMsMSxmdW5jdGlvbih4LmkpIGMoeC5pICUqJSBBICUqJSB4LmkpICsgc3VtKHguaSAqICgyKnEvc3VtLnEgKyBjKEEgJSolIG9tZWdhMikpKSkgKyBwbm9ybSh0YXUudGlsZGUsbG9nLnA9VFJVRSkKICBpZihsb2cpCiAgICByZXR1cm4odmFsKQogIGVsc2UKICAgIHJldHVybihleHAodmFsKSkgICAgCn0KCmRsb2dfYmV0YSAgPC0gZnVuY3Rpb24oeCwgYSwgYikgZGJldGEoeCwgYSwgYiwgbG9nID0gVFJVRSkgIyBtaXh0dXJlIHdlaWdodCBwcmlvcgpkbG9nX2dhbW1hIDwtIGZ1bmN0aW9uKHgsIGEsIGIpIGRnYW1tYSh4LCBzaGFwZSA9IGEsIHJhdGUgPSBiLCBsb2cgPSBUUlVFKQoKY29tcG9uZW50X2xvZ2xpayA8LSBmdW5jdGlvbihYLCBjb29yZCwgYWxwaGEsIHRoZXRhLCBhbmNob3IgPSAxKXsKICBHYW1tYSA8LSBidWlsZF9HYW1tYShjb29yZCA9IGNvb3JkLCBhbHBoYSA9IGFscGhhLCB0aGV0YSA9IHRoZXRhKQogIHByZXAgIDwtIHByZXBfYnIoR2FtbWEpCiAgaWR4ICAgPC0gcHJlcFtbYW5jaG9yXV0kaWR4CiAgTCAgICAgPC0gcHJlcFtbYW5jaG9yXV0kTAogIFNpZ21hIDwtIHQoTCkgJSolIEwKICBYc3ViICA8LSBhcy5tYXRyaXgoWFssIGlkeCwgZHJvcCA9IEZBTFNFXSkKICBpbnRlbnNpdHlfbG9nc2tldygKICAgIFhzdWIsIHBhciA9IGxpc3QoU2lnbWEsIGRlbHRhID0gcmVwKDAsIG5jb2woWHN1YikpKSwKICAgIGFscGhhLnBhcmEgPSBGQUxTRSwgbG9nID0gVFJVRQogICkKfQoKbWl4X2xvZ2xpayA8LSBmdW5jdGlvbihYLCBjb29yZCwgYWxwaGEsIHRoMSwgdGgyLCBhbmNob3IgPSAxKSB7CiAgaWYgKCFpcy5maW5pdGUodGgxKSB8fCB0aDEgPD0gMCkgcmV0dXJuKGxpc3QobDEgPSByZXAoLUluZiwgbnJvdyhYKSksIGwyID0gcmVwKC1JbmYsIG5yb3coWCkpKSkKICBpZiAoIWlzLmZpbml0ZSh0aDIpIHx8IHRoMiA8PSAwKSByZXR1cm4obGlzdChsMSA9IHJlcCgtSW5mLCBucm93KFgpKSwgbDIgPSByZXAoLUluZiwgbnJvdyhYKSkpKQogIGwxIDwtIGNvbXBvbmVudF9sb2dsaWsoWCwgY29vcmQsIGFscGhhID0gYWxwaGEsIHRoZXRhID0gdGgxLCBhbmNob3IgPSBhbmNob3IpCiAgbDIgPC0gY29tcG9uZW50X2xvZ2xpayhYLCBjb29yZCwgYWxwaGEgPSBhbHBoYSwgdGhldGEgPSB0aDIsIGFuY2hvciA9IGFuY2hvcikKICBsaXN0KGwxID0gbDEsIGwyID0gbDIpCn0KCnRhdV9mcm9tX2xsIDwtIGZ1bmN0aW9uKGwxLCBsMiwgbGFtYmRhKSB7CiAgYSA8LSBsb2cobGFtYmRhKSAgICAgKyBsMQogIGIgPC0gbG9nKDEgLSBsYW1iZGEpICsgbDIKICBtIDwtIHBtYXgoYSwgYikKICBleHAoYSAtIG0pIC8gKGV4cChhIC0gbSkgKyBleHAoYiAtIG0pKQp9CgptY21jX2JyIDwtIGZ1bmN0aW9uKAogICAgZGF0YSwgY29vcmQsIGFscGhhLAogICAgbl9pdGVyICAgPSA0MDAwLAogICAgYnVybiAgICAgPSAxMDAwLAogICAgdGhpbiAgICAgPSAyLAogICAgIyBwcmlvcnMKICAgIGFfbGFtYmRhID0gMSwgYl9sYW1iZGEgPSAxLCAgICMgZm9yIG1peHR1cmUgd2VpZ2h0IGxhbWJkYQogICAgYV90aDIgPSAxLCBiX3RoMiA9IDEsICAgICAgICAgIyBmb3IgdGhldGEyCiAgICBhX2RlbCA9IDEsIGJfZGVsID0gMSwgICAgICAgICAjIGZvciBkZWx0YQogICAgIyBwcm9wb3NhbCBzZHMgb24gbG9nLXNjYWxlCiAgICBwcm9wX3NkX2xvZ190aDIgID0gMC4xMCwKICAgIHByb3Bfc2RfbG9nX2RlbCAgPSAwLjEwLAogICAgIyBpbml0aWFscwogICAgaW5pdF9sYW1iZGEgPSAwLjUsICAjIG1peHR1cmUgd2VpZ2h0CiAgICBpbml0X3RoMiAgICA9IDE1LCAgICMgdGhldGEyCiAgICBpbml0X2RlbCAgICA9IDIwLCAgICMgZGVsdGEKICAgIGFuY2hvciA9IDEsCiAgICB2ZXJib3NlX2V2ZXJ5ID0gMTAwCil7CiAgc3RvcGlmbm90KGluaXRfbGFtYmRhID4gMCwgaW5pdF9sYW1iZGEgPCAxLCBpbml0X3RoMiA+IDAsIGluaXRfZGVsID4gMCkKICAKICBrZWVwX2lkeCA8LSBzZXEuaW50KGJ1cm4gKyAxLCBuX2l0ZXIsIGJ5ID0gdGhpbikKICBuX2tlZXAgICA8LSBsZW5ndGgoa2VlcF9pZHgpCiAgCiAgZHJhd3MgPC0gZGF0YS5mcmFtZSgKICAgIGl0ZXIgICAgID0ga2VlcF9pZHgsCiAgICBsYW1iZGEgICA9IE5BX3JlYWxfLCAgICMgbWl4dHVyZSB3ZWlnaHQKICAgIHRoZXRhMSAgID0gTkFfcmVhbF8sCiAgICB0aGV0YTIgICA9IE5BX3JlYWxfCiAgKQogIHpfZHJhd3MgPC0gbWF0cml4KE5BX2ludGVnZXJfLCBucm93ID0gbl9rZWVwLCBuY29sID0gbnJvdyhkYXRhKSkKICB0YXVfYmFyIDwtIHJlcCgwLCBucm93KGRhdGEpKQogIAogICMgaW5pdAogIGxhbWJkYSA8LSBpbml0X2xhbWJkYQogIHRoMiAgICA8LSBpbml0X3RoMgogIGRlbCAgICA8LSBpbml0X2RlbAogIHRoMSAgICA8LSB0aDIgKyBkZWwKICAKICBsbCAgPC0gbWl4X2xvZ2xpayhYID0gZGF0YSwgY29vcmQgPSBjb29yZCwgYWxwaGEgPSBhbHBoYSwgdGgxID0gdGgxLCB0aDIgPSB0aDIsIGFuY2hvciA9IGFuY2hvcikKICB0YXUgPC0gdGF1X2Zyb21fbGwobGwkbDEsIGxsJGwyLCBsYW1iZGEgPSBsYW1iZGEpCiAgeiAgIDwtIHJiaW5vbShuID0gbGVuZ3RoKHRhdSksIHNpemUgPSAxLCBwcm9iID0gdGF1KQogIAogIGFjY19ibG9jayA8LSAwTDsgbl9wcm9wIDwtIDBMCiAgCiAgIyBoZWxwZXI6IGNvbXBsZXRlLWRhdGEgbG9nLXBvc3RlcmlvciBpbiDOty1zcGFjZSAozrcyPWxvZyDOuDIsIM63ZD1sb2cgzrQpCiAgbG9ncG9zdF9ldGEgPC0gZnVuY3Rpb24oZXRhMiwgZXRhZCwgel92ZWMpIHsKICAgIHRoMnAgPC0gZXhwKGV0YTIpCiAgICBkZWxwIDwtIGV4cChldGFkKQogICAgdGgxcCA8LSB0aDJwICsgZGVscAogICAgCiAgICBsbHAgPC0gbWl4X2xvZ2xpayhYID0gZGF0YSwgY29vcmQgPSBjb29yZCwgYWxwaGEgPSBhbHBoYSwKICAgICAgICAgICAgICAgICAgICAgIHRoMSA9IHRoMXAsIHRoMiA9IHRoMnAsIGFuY2hvciA9IGFuY2hvcikKICAgIAogICAgbHBfbGwgICAgPC0gc3VtKHpfdmVjICogbGxwJGwxICsgKDEgLSB6X3ZlYykgKiBsbHAkbDIpCiAgICBscF9wcmlvciA8LSBkbG9nX2dhbW1hKHRoMnAsIGFfdGgyLCBiX3RoMikgKyBkbG9nX2dhbW1hKGRlbHAsIGFfZGVsLCBiX2RlbCkKICAgIGxwX2phYyAgIDwtIGV0YTIgKyBldGFkCiAgICBsaXN0KHZhbHVlID0gbHBfbGwgKyBscF9wcmlvciArIGxwX2phYywKICAgICAgICAgdGgxID0gdGgxcCwgdGgyID0gdGgycCwgZGVsID0gZGVscCwgbDFfdmVjID0gbGxwJGwxLCBsMl92ZWMgPSBsbHAkbDIpCiAgfQogIAogIGtlZXBfY3Vyc29yIDwtIDFMCiAgZm9yIChpdCBpbiBzZXFfbGVuKG5faXRlcikpIHsKICAgIAogICAgIyMgKDEpIGJsb2NrZWQgTUggZm9yICjOuDIsIM60KSBvbiBsb2ctc2NhbGUKICAgIG5fcHJvcCA8LSBuX3Byb3AgKyAxTAogICAgZXRhMiAgPC0gbG9nKHRoMik7ICBldGFkIDwtIGxvZyhkZWwpCiAgICBldGEycCA8LSBybm9ybSgxLCBtZWFuID0gZXRhMiwgc2QgPSBwcm9wX3NkX2xvZ190aDIpCiAgICBldGFkcCA8LSBybm9ybSgxLCBtZWFuID0gZXRhZCwgc2QgPSBwcm9wX3NkX2xvZ19kZWwpCiAgICAKICAgICMgY3VycmVudCB2YWx1ZSB2aWEgc2FtZSBmdW5jdGlvbgogICAgY3VyciA8LSBsb2dwb3N0X2V0YShldGEyLCBldGFkLCB6KQogICAgcHJvcCA8LSBsb2dwb3N0X2V0YShldGEycCwgZXRhZHAsIHopCiAgICAKICAgIGlmIChsb2cocnVuaWYoMSkpIDwgKHByb3AkdmFsdWUgLSBjdXJyJHZhbHVlKSkgewogICAgICB0aDIgPC0gcHJvcCR0aDI7ICBkZWwgPC0gcHJvcCRkZWw7ICB0aDEgPC0gcHJvcCR0aDEKICAgICAgbGwkbDEgPC0gcHJvcCRsMV92ZWM7IGxsJGwyIDwtIHByb3AkbDJfdmVjCiAgICAgIGFjY19ibG9jayA8LSBhY2NfYmxvY2sgKyAxTAogICAgfQogICAgIyBlbHNlIGtlZXAgKHRoMSx0aDIsZGVsLGxsKSBhcy1pcwogICAgCiAgICAjIyAoMikgR2liYnMgZm9yIGxhbWJkYSB8IHogKEJldGEtY29uanVnYXRlKQogICAgbjEgPC0gc3VtKHopOyBuMCA8LSBsZW5ndGgoeikgLSBuMQogICAgbGFtYmRhIDwtIHJiZXRhKDEsIGFfbGFtYmRhICsgbjEsIGJfbGFtYmRhICsgbjApCiAgICAKICAgICMjICgzKSBHaWJicyBmb3IgeiB8IHgsIM64LCDOuwogICAgdGF1IDwtIHRhdV9mcm9tX2xsKGxsJGwxLCBsbCRsMiwgbGFtYmRhKQogICAgeiAgIDwtIHJiaW5vbShsZW5ndGgodGF1KSwgMSwgdGF1KQogICAgCiAgICAjIyBzdG9yZSB0aGlubmVkIGRyYXdzCiAgICBpZiAoaXQgJWluJSBrZWVwX2lkeCkgewogICAgICB0YXVfYmFyIDwtIHRhdV9iYXIgKyB0YXUgLyBuX2tlZXAKICAgICAgZHJhd3MkbGFtYmRhW2tlZXBfY3Vyc29yXSA8LSBsYW1iZGEKICAgICAgZHJhd3MkdGhldGExW2tlZXBfY3Vyc29yXSA8LSB0aDEKICAgICAgZHJhd3MkdGhldGEyW2tlZXBfY3Vyc29yXSA8LSB0aDIKICAgICAgel9kcmF3c1trZWVwX2N1cnNvciwgXSAgICA8LSB6CiAgICAgIGtlZXBfY3Vyc29yIDwtIGtlZXBfY3Vyc29yICsgMUwKICAgIH0KICAgIAogICAgaWYgKHZlcmJvc2VfZXZlcnkgPiAwICYmIGl0ICUlIHZlcmJvc2VfZXZlcnkgPT0gMCkgewogICAgICBjYXQoc3ByaW50ZigiaXRlciAlZCB8IGxhbWJkYT0lLjNmIM64MT0lLjNmIM64Mj0lLjNmIHwgYWNjKM64LWJsb2NrKT0lLjNmXG4iLAogICAgICAgICAgICAgICAgICBpdCwgbGFtYmRhLCB0aDEsIHRoMiwgYWNjX2Jsb2NrIC8gbl9wcm9wKSkKICAgIH0KICB9CiAgCiAgbGlzdCgKICAgIGRyYXdzID0gZHJhd3MsCiAgICB6X2RyYXdzID0gel9kcmF3cywKICAgIGFjY2VwdF9yYXRlX3RoZXRhX2Jsb2NrID0gYWNjX2Jsb2NrIC8gbl9wcm9wLAogICAgdGF1X2JhciA9IHRhdV9iYXIsCiAgICBoYXJkX2xhYmVscyA9IGlmZWxzZSh0YXVfYmFyID4gMC41LCAxLCAyKSwKICAgIHBvc3RfbWVhbiA9IGMobGFtYmRhID0gbWVhbihkcmF3cyRsYW1iZGEpLAogICAgICAgICAgICAgICAgICB0aGV0YTEgPSBtZWFuKGRyYXdzJHRoZXRhMSksCiAgICAgICAgICAgICAgICAgIHRoZXRhMiA9IG1lYW4oZHJhd3MkdGhldGEyKSksCiAgICBwb3N0X21lZGlhbiA9IGMobGFtYmRhID0gbWVkaWFuKGRyYXdzJGxhbWJkYSksCiAgICAgICAgICAgICAgICAgICAgdGhldGExID0gbWVkaWFuKGRyYXdzJHRoZXRhMSksCiAgICAgICAgICAgICAgICAgICAgdGhldGEyID0gbWVkaWFuKGRyYXdzJHRoZXRhMikpCiAgKQp9CmBgYAoKZml0IHRoZSBkYXRhCmBgYHtyfQphbHBoYSA8LSAxLjUKCmZpdDEgPC0gbWNtY19icigKICBkYXRhID0gZXh0cmVtZV9MMV9ob3VybHksIGNvb3JkID0gY29vcmQsIGFscGhhID0gYWxwaGEsCiAgbl9pdGVyID0gMTAwMDAsIGJ1cm4gPSAyMDAwLCB0aGluID0gMiwKICAjIHByaW9ycwogIGFfbGFtYmRhID0gMSwgYl9sYW1iZGEgPSAxLCAgICMgcHJpb3IgZm9yIG1peHR1cmUgd2VpZ2h0IGxhbWJkYQogIGFfdGgyID0gMSwgYl90aDIgPSAxLCAgICAgICAgICMgcHJpb3IgZm9yIHRoZXRhMgogIGFfZGVsID0gMSwgYl9kZWwgPSAxLCAgICAgICAgICMgcHJpb3IgZm9yIGRlbHRhCiAgIyBpbml0aWFscwogIGluaXRfdGgyICA9IDIsCiAgaW5pdF9kZWwgID0gNCwKICAjIHByb3Bvc2FscyAobG9nLXNjYWxlKQogIHByb3Bfc2RfbG9nX3RoMiAgPSAwLjAxLAogIHByb3Bfc2RfbG9nX2RlbCAgPSAwLjAxLAogIGFuY2hvciA9IDEsCiAgdmVyYm9zZV9ldmVyeSA9IDUwCikKCnNhdmUoZml0MSwgZmlsZSA9ICJhcHBsaWNhdGlvbl9maXRfaG91cmx5X2ZpbmFsLlJEYXRhIikKYGBgCgpgYGB7cn0KZGYgPC0gZml0MSRkcmF3cwpnMSA8LSBnZ3Bsb3QoZGYsIGFlcyhpdGVyLCB0aGV0YTEpKSArIGdlb21fbGluZSgpICsKICBsYWJzKHRpdGxlID0gZXhwcmVzc2lvbih0cmFjZX50aGV0YVsxXSkpCmcyIDwtIGdncGxvdChkZiwgYWVzKGl0ZXIsIHRoZXRhMikpICsgZ2VvbV9saW5lKCkgKwogIGxhYnModGl0bGUgPSBleHByZXNzaW9uKHRyYWNlfnRoZXRhWzJdKSkKZzMgPC0gZ2dwbG90KGRmLCBhZXMoaXRlciwgbGFtYmRhKSkgKyBnZW9tX2xpbmUoKSArCiAgbGFicyh0aXRsZSA9IGV4cHJlc3Npb24odHJhY2V+bGFtYmRhKSkKZ2QxIDwtIGdncGxvdChkZiwgYWVzKHRoZXRhMSkpICsgZ2VvbV9kZW5zaXR5KCkgKwogIGxhYnModGl0bGUgPSBleHByZXNzaW9uKGRlbnNpdHl+dGhldGFbMV0pKQpnZDIgPC0gZ2dwbG90KGRmLCBhZXModGhldGEyKSkgKyBnZW9tX2RlbnNpdHkoKSArCiAgbGFicyh0aXRsZSA9IGV4cHJlc3Npb24oZGVuc2l0eX50aGV0YVsyXSkpCmdkMyA8LSBnZ3Bsb3QoZGYsIGFlcyhsYW1iZGEpKSArIGdlb21fZGVuc2l0eSgpICsKICBsYWJzKHRpdGxlID0gZXhwcmVzc2lvbihkZW5zaXR5fmxhbWJkYSkpCgpncmlkRXh0cmE6OmdyaWQuYXJyYW5nZShnMSxnMixnMywgZ2QxLGdkMixnZDMsIG5jb2wgPSAzKQpgYGAKYGBge3J9CmZpdDEkcG9zdF9tZWFuCmBgYAoKYGBge3J9CnNkKGZpdDEkZHJhd3MkdGhldGExKQpzZChmaXQxJGRyYXdzJHRoZXRhMikKYGBgCgoKYGBge3J9CnBsb3RfZXh0cmVtZV9vbl9tYXAgPC0gZnVuY3Rpb24oaSwgZXh0cmVtZV9tYXQgPSBleHRyZW1lX0wxX2hvdXJseSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtYXBfc2YgPSBzcXVhcmVfbWFwLCBkYXRhX3dpZGUgPSBkYXRhX2hvdXJseV9zdWIpIHsKICBzdG9waWZub3QoaXMubWF0cml4KGV4dHJlbWVfbWF0KSwKICAgICAgICAgICAgaSA+PSAxLCBpIDw9IG5yb3coZXh0cmVtZV9tYXQpKQogIAogICMgY29sdW1uIG9yZGVyIG9mIHBpeGVscyBpbiB0aGUgd2lkZSBkYWlseSBkYXRhIHVzZWQgdG8gYnVpbGQgcmFpbl9tYXQvZXh0cmVtZV9MMQogIHBpeF9vcmRlciA8LSBhcy5pbnRlZ2VyKHN1YigiUCIsICIiLCBjb2xuYW1lcyhkYXRhX3dpZGUpWy0xXSkpCiAgCiAgIyB2YWx1ZXMgZm9yIHRoZSBpLXRoIGV4dHJlbWUgZXZlbnQKICB2YWxzX2kgPC0gYXMubnVtZXJpYyhleHRyZW1lX21hdFtpLCBdKQogIGRmX2kgPC0gZGF0YS5mcmFtZShQSVhFTCA9IHBpeF9vcmRlciwgdmFsX2kgPSB2YWxzX2kpCiAgCiAgIyBqb2luIHRvIGdlb21ldHJ5CiAgam9pbmVkIDwtIG1hcF9zZiAlPiUKICAgIGxlZnRfam9pbihkZl9pLCBieSA9ICJQSVhFTCIpCiAgCiAgIyBxdWljayBzdHJldGNoIGZvciBuaWNlciBsZWdlbmQgYnJlYWtzCiAgYnJrcyA8LSBxdWFudGlsZShkZl9pJHZhbF9pLCBwcm9icyA9IHNlcSgwLCAxLCBsZW5ndGgub3V0ID0gNiksIG5hLnJtID0gVFJVRSkKICBicmtzX2ZpeGVkIDwtIHNvcnQodW5pcXVlKGJya3MpKQogICMgcGxvdAogIG1hcHZpZXcoCiAgICBqb2luZWQsCiAgICB6Y29sICAgICAgID0gInZhbF9pIiwKICAgIGxheWVyLm5hbWUgPSBwYXN0ZTAoIkV4dHJlbWUgTDEgc25hcHNob3QgKHJvdyAiLCBpLCAiKSIpLAogICAgYXQgICAgICAgICA9IGJya3NfZml4ZWQKICApCn0KYGBgCgpgYGB7cn0KcGxvdF9leHRyZW1lX29uX21hcCgxKQpgYGAKCmBgYHtyfQpwbG90X2V4dHJlbWVfb25fbWFwKDUpCmBgYAoKJFxjaGleMihTX0EsIFNfQikkIHRlc3QKYGBge3J9CmNoaV9xIDwtIGZ1bmN0aW9uKFgsIHEgPSAwLjk4KSB7CiAgbiA8LSBuY29sKFgpCiAgbSA8LSBucm93KFgpCiAgCiAgIyB0aHJlc2hvbGRzIGF0IHNpdGUgbGV2ZWwKICB1IDwtIGFwcGx5KFgsIDIsIHF1YW50aWxlLCBwcm9icyA9IHEsIG5hLnJtID0gVFJVRSkKICAKICAjIGluaXRpYWxpemUgbWF0cml4CiAgY2hpX2hhdCA8LSBtYXRyaXgoTkEsIG4sIG4pCiAgCiAgZm9yIChpIGluIDE6KG4tMSkpIHsKICAgIGZvciAoaiBpbiAoaSsxKTpuKSB7CiAgICAgIEFpIDwtIFhbLCBpXSA+IHVbaV0KICAgICAgQWogPC0gWFssIGpdID4gdVtqXQogICAgICBjaGlfaGF0W2ksIGpdIDwtIG1lYW4oQWkgJiBBaikgLyBtZWFuKEFpKQogICAgICBjaGlfaGF0W2osIGldIDwtIGNoaV9oYXRbaSwgal0KICAgIH0KICB9CiAgZGlhZyhjaGlfaGF0KSA8LSAxCiAgY2hpX2hhdAp9CgpkaXN0X21hdHJpeCA8LSBhcy5tYXRyaXgoZGlzdChjb29yZCkpCgpzdG9waWZub3QobGVuZ3RoKGZpdDEkaGFyZF9sYWJlbHMpID09IG5yb3coZXh0cmVtZV9MMV9ob3VybHkpKQpsYWJzIDwtIGZpdDEkaGFyZF9sYWJlbHMKClgxIDwtIGV4dHJlbWVfTDFfaG91cmx5W2xhYnMgPT0gMSwgLCBkcm9wID0gRkFMU0VdClgyIDwtIGV4dHJlbWVfTDFfaG91cmx5W2xhYnMgPT0gMiwgLCBkcm9wID0gRkFMU0VdCgojIDQpIGNvbXB1dGUgz4fMgl9xIGZvciBib3RoIGNvbXBzIGFuZCB0d28gcSBsZXZlbHMKY2hpXzk1X2MxIDwtIGNoaV9xKFgxLCBxID0gMC45NSkKY2hpXzk1X2MyIDwtIGNoaV9xKFgyLCBxID0gMC45NSkKY2hpXzkwX2MxIDwtIGNoaV9xKFgxLCBxID0gMC45MCkKY2hpXzkwX2MyIDwtIGNoaV9xKFgyLCBxID0gMC45MCkKYGBgCgoKYGBge3J9CnVwcGVyX3RvX2RmIDwtIGZ1bmN0aW9uKGNoaV9tYXQsIGRpc3RfbWF0LCBxX2xhYmVsLCBjb21wX2xhYmVsKSB7CiAgaWR4IDwtIHdoaWNoKHVwcGVyLnRyaShjaGlfbWF0KSwgYXJyLmluZCA9IFRSVUUpCiAgdGliYmxlKAogICAgZGlzdCA9IGRpc3RfbWF0W2lkeF0sCiAgICBjaGkgID0gY2hpX21hdFtpZHhdLAogICAgcSAgICA9IHFfbGFiZWwsCiAgICBjb21wID0gY29tcF9sYWJlbAogICkgJT4lIGZpbHRlcihpcy5maW5pdGUoY2hpKSkKfQoKZGZfYWxsIDwtIGJpbmRfcm93cygKICB1cHBlcl90b19kZihjaGlfOTVfYzEsIGRpc3RfbWF0cml4LCAicT0wLjk1IiwgIkNvbXAxIiksCiAgdXBwZXJfdG9fZGYoY2hpXzk1X2MyLCBkaXN0X21hdHJpeCwgInE9MC45NSIsICJDb21wMiIpLAogIHVwcGVyX3RvX2RmKGNoaV85MF9jMSwgZGlzdF9tYXRyaXgsICJxPTAuOTAiLCAiQ29tcDEiKSwKICB1cHBlcl90b19kZihjaGlfOTBfYzIsIGRpc3RfbWF0cml4LCAicT0wLjkwIiwgIkNvbXAyIiksCgopCgojIDYpIHBsb3QgKG1pcnJvcnMgeW91ciBvcmlnaW5hbCBzdHlsZSkKZ2dwbG90KGRmX2FsbCwgYWVzKHggPSBkaXN0LCB5ID0gY2hpLCBjb2xvciA9IGNvbXApKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuNCwgc2l6ZSA9IDAuNikgKwogIGZhY2V0X3dyYXAofiBxLCBzY2FsZXMgPSAiZnJlZV95IikgKwogIHRoZW1lX2J3KCkgKwogIGxhYnMoCiAgICB4ID0gZXhwcmVzc2lvbihoKHNbQV0sIHNbQl0pKSwKICAgIHkgPSBleHByZXNzaW9uKGhhdChjaGkpW3FdKHNbQV0sIHNbQl0pKSwKICAgIHRpdGxlID0gIkVtcGlyaWNhbCAiICUrJSBleHByZXNzaW9uKGhhdChjaGkpW3FdKSAlKyUgIiB2cyBEaXN0YW5jZSBieSBNQ01DIEhhcmQtTGFiZWwgQ29tcG9uZW50IgogICkgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCJDb21wMSIgPSAiYmx1ZSIsICJDb21wMiIgPSAicmVkIikpCgpgYGAKCgoKYGBge3J9CmZvciAoaSBpbiBzZXFfbGVuKG5jb2woZml0MSR6X2RyYXdzKSkpIHsKICBwbG90KAogICAgZml0MSR6X2RyYXdzWywgaV0sCiAgICB0eXBlID0gImwiLAogICAgeGxhYiA9ICJpdGVyYXRpb25zIiwKICAgIHlsYWIgPSAieiBkcmF3IiwKICAgIG1haW4gPSBwYXN0ZSgic2FtcGxlIiwgaSkgIyA8LS0tIE1vZGlmaWVkIGxpbmUKICApCn0KYGBgCmBgYHtyfQoKCnBsb3QodGltZXNfaG91cmx5LCBmaXQxJHRhdV9iYXIsIHR5cGUgPSAibCIsCiAgICAgeGxhYiA9ICJUaW1lIiwgeWxhYiA9ICJQb3N0ZXJpb3IgbWVhbiBvZiB6IiwKICAgICBtYWluID0gIlBvc3RlcmlvciBtZWFuIHdpdGggc2Vhc29uIikKCmBgYAoKCgoK